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1 .  Introduction 

---In  general,  the  objective  of  a  computational  solution  of  a  practical  prob¬ 
lem  is  to  obtain  an  acceptably  accurate  and  reliable  prediction  of  the  behavior 
of  the  physical  phenomena  under  study.  For  this  the  computation  should  produce 
not  only  an  approximate  numerical  solution  but  also  some  reliable  information 
about  its  accuracy.  , 

The  selection  of  the  specific  accuracy  requirement  depends  strongly  on  the 
goal  of  the  computation.  Often  it  is  desired  to  obtain  detailed  information 
about  the  solution  itself.  In  other  cases,  the  main  focus  is  the  value  of  a 
specified  functional  of  the  solution,  as,  for  instance,  a  stress  intensity  fac¬ 
tor  in  fracture  mechanics,  or  a  drag  in  fluid  dynamics.  Then,  again  in  other 
cases,  interest  centers,  say,  on  streamlines  and/or  on  the  existence  of  separa¬ 
tion  lines,  as,  for  example,  in  sepage  problems.  Other  goals  may  be  the  deter¬ 
mination  of  certain  critical  data,  as,  collapse  points  or  buckling  points  and 
their  associated  loads  in  structural  mechanics,  etc. 

In  connection  with  most  of  these  goals  we  are  interested  in  quantitative 
results  which  have  a  desired  accuracy.  For  this  the  error  has  to  be  defined, 
that  is,  the  exact  results  have  to  be  specified  with  which  the  computed  data  are 
to  be  compared,  a  norm  has  to  be  prescribed  under  which  the  error  is  to  be  mea¬ 
sured,  and,  last  but  not  least,  some  procedure  has  to  be  established  for  estim- 

1)  This  work  was  in  part  supported  by  ONR-contracts  N  0014-77-C-0623 
N  0014-80-C-04  55 . 
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lege  Park,  MD  20742. 

3)  Department  of  Math  and  Statistics,  Univ.  of  Pittsburgh,  Pittsburgh,  PA  15238. 
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ating  the  desired  error-data.  On  the  other  hand,  for  some  of  the  indicated  goals 
we  may  he  interested  only  in  qualitative  behavior,  as  for  instance  the  mentioned 
existence  of  separation  lines.  Such  a  requirement  may  lead  into  other  branches 
of  mathematics  as,  for  example,  into  the  qualitative  theory  of  ordinary  differ¬ 
ential  equations  in  the  case  of  streamline  problems. 

In  recent  years,  the  use  of  the  finite-element  in  engineering  and  science 
has  become  a  multi-million  dollar  enterprise.  Many  programs  exist  and  many  more 
are  being  developed,  (see  eq.  [1]).  In  line  with  the  stated  need  for  information 
about  the  accuracy  of  the  computed  solution  we  may  expect  that  future  programs 
will  incorporate  facilities  for  computing  such  error  estimates.  As  indicated 
before,  these  estimates  should  correspond  to  the  principal  goals  of  the  computa¬ 
tion  that  can  be  selected  by  the  user  of  the  program,  and  the  program  should 
allow  a  choice  of  the  norm  under  which  the  estimate  is  to  be  computed. 

Such  error-estimation  facilities  are  certainly  of  direct  importance  in  many 
applications,  provided  they  can  be  guaranteed  to  be  reliable.  For  instance,  for 
the  many  types  of  certif ications-computations  required  in  the  design  of  complex 
structures,  nuclear  plants  etc.,  the  availability  of  reliable  estimates  of  the 
accuracy  of  the  computed  data  is  obviously  very  essential.  In  other  cases,  such 
estimates  may  reduce  the  total  design  effort  and  avoid  unnecessary  overdesign. 

On  the  other  hand,  the  availability  of  effective  error  estimates  introduces 
the  possibility  of  structuring  the  finite  element  computation  so  as  to  achieve 
a  desired  error  tolerance  at  minimal  cost  or  to  provide  the  best  possible  solu¬ 
tion  within  an  allowable  cost  range.  This  is  an  important  objective  especially 
in  view  of  the  fact  that  in  typical  engineering  applications  an  accuracy  of,  say, 
2-5%  is  desirable,  and  any  effort  expended  to  obtain  a  much  higher  accuracy  is 
a  waste. 

It  is  beginning  to  be  a  widely  accepted  observation  that  for  realistic 
problems  it  is  rarely  feasible  to  design  numerical  processes  which  reliably  and 
effectively  achieve  a  desired  accuracy  at  reasonable  cost  and  yet  which  do  not 
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utilize  some  form  of  adaptivity.  At  the  same  time  adaptive  approaches  may  also 
lead  to  a  simplification  of  the  input  data  needed  for  the  program  and  hence  free 
the  user  of  part  of  t he  drudgery  typical  in  preparing  such  input  and  of  relieving 
ines  cf  manv  of  the  a  priori  decisions  required  by  most  of  today's  programs. 

While  the  need  for  adaptive  aoproaches  is  well  recognized,  the  precise 
meaning  of  the  adaptation-concept  remains  vague  and  ill-defined.  A  classifica¬ 
tion  of  this  concept  may  require  close  attention  to  results  in  such  fields  as 
automatic  control  theory,  artificial  intelligence,  learning  processes,  etc.  (see 
eq.  [2],  [3],  [4],  [5]).  In  our  context,  the  availability  of  precisely  defined 
reliable  error  estmators  appears  to  be  central  to  the  design  of  effective  adap¬ 
tive  procedures.  The  effectivity  has  to  be  measured  here  in  terms  of  the  size 
of  the  iass  of  prot lems  for  which  solutions  within  a  specified  range  of  accuracy 
are  obtainable  with  minimal  cost.  Of  course,  the  question  of  cost  is  here  a  very 
complex  issue  and  still  require?  much  thought. 

By  far  tne  majority  of  todays  finite  element  computations  involve  linear 
problems,  altnough  more  and  more  attention  is  being  paid  to  nonlinear  phenomena. 
Not  unexpectedly,  there  are  many  differences  between  linear  and  nonlinear  pro¬ 
blems,  and  the  latter  show  many  special  features  not  present  in  the  linear  case. 
\s  stated  at  the  outset,  the  objective  of  the  computation  should  be  a  prediction 
of  the  behavior  of  the  physical  system  under  study.  If  we  are  concerned, for 
instance,  with  a  mechanical  structure,  then  the  term  "behavior"  actually  refers, 
say,  to  the  deformations  of  the  structure  in  response  to  the  action  of  external 
influence  quantities,  such  as  loads,  changes  in  material  properties  and  geometri¬ 
cal  data,  etc.  This  means  that  the  equations  which  govern  the  deformation  o  ■ 
behavior  variables  also  depend  on  a  set  of  influence  parameters  and  the  objective 
is  to  assess  the  behavior  of  the  solutions  under  general  variations  of  the  para¬ 
meters.  In  the  linear  case,  most  of  this  variation  is  also  assumed  to  be  linear 
and  we  need  only  compute  a  few  specific  solutions  to  achieve  our  objective. 


In  nonlineal  problems,  however,  the  computation  of  a  few  specific  solutions  pro¬ 
vides  little  or  no  insight  into  the  behavior  of  the  system  and  we  are  led  to 
onsider  the  set  of  all  solutions  depending  on  all  the  parameters  in  a  specific 
range.  In  structural  mechanics,  this  set  is  now  usually  called  the  equilibrium 
surface  of  the  structure.  The  objective  then  becomes  an  analysis  of  the  form 
and  special  features  of  this  surface,  as  for  example,  the  boundaries  of  the 
stability  regions  on  the  surface,  etc.  The  principal  tools  for  such  an  analysis 
are  the  continuation  methods,  as  incremental  processes,  as  they  are  called  in 
structural  mechanics.  As  before  error  estimations  and  adaptive  approaches  will 
have  to  be  incorporated  into  these  processes.  This  is  as  yet  a  wide-open  re¬ 
search  area. 

This  survey  is  intended  to  give  an  overview  of  some  of  the  recent  results 
that  have  become  available  about  reliable,  computable  error  estimators  for  finite 
element  solutions  as  they  exist  in  the  design  of  adaptive  solution  procedures 
both  in  the  linear  as  well  as  the  nonlinear  case.  By  necessity  the  discussion 
had  to  be  kept  on  a  descriptive  level  and  for  all  mathematical  details  we  need 
to  refer  to  the  cited  literature. 


2.  Accuracy  Considerations  in  the  Finite  Clement  Method 


Tiie  finite  element  approach  to  the  solutions  of  suitable  boundary  value 
problems  raav  be  cl aracterized  loosely  hv  the  use  of  a  weak  bivariate  form  of 
the  problem  and  the  approximation  of  the  solution  by  piecewise  smooth  functions. 
There  is  a  growing  number  of  introductions  to  the  basic  theory  of  the  method; 
we  refer  only  to  [6],  [7],  [8],  and  [9]  to  mention  a  few. 

Many  of  the  available  theoretical  results  concern  an  asymptotic  analysis  of 
the  convergence  rate  of  the  method.  Although  this  type  of  analys  s  is  not  too 
useful  for  assessing  the  accuracy  of  a  specific  solution,  it  provides  important 
insight  and  offers  guidelines  for  future  efforts. 

We  illustrate  this  point  with  some  observations  from  general  practical  exper¬ 
ience  with  engineering  problems.  It  is  well  accepted  that  in  the  case  of  fairly 
complex  problems  a  relative  accuracy  of  about  8-10%  in  the  energy  norm  is  achievable 
with  about  1000  to  2000  degrees  of  freedom  N.  On  the  other  hand,  in  order  to  obtain 
an  accuracy  of  3-5%,  the  rate  of  convergence  should  be  at  least  of  the  order  of 
N  But.  in  the  presence  of  singularities  and  for  uniform  meshes  in  two  space 
dimensions  the  "standard  method"  provides  often  only  a  rate  of  0(1,  *  or  worse,  and 
hence  the  desired  accuracy  of  3-5%  is  practically  infeasible.  For  three-dimensional 
problems  the  situation  is  even  more  difficult. 

By  the  standard  method  we  mean  here  the  so-called  h-version  of  the  finite 
element  method  which  is  most  widely  used  and  studied.  In  that  version,  the  degrees 
of  the  elements  used  for  the  computation  are  fixed  and  the  increasing  accuracy  is 
achieved  by  the  use  of  appropriate  sequences  of  finer  and  finer  meshes.  The  In¬ 
version  indeed  underlies  the  design  of  practically  all  of  today's  commercially 
available  codes  (see  [1]). 

Recently  a  new  approach  became  available  with  the  development  of  the  so-called 
p-version  and  h-p-version  of  the  finite  element  method.  In  the  p-version,  the 
meshes  are  fixed  and  the  increasing  accuracy  is  achieved  by  the  use  of  elements 
with  higher  and  higher  degrees.  The  h-p-version  combines  the  two  approaches  in  a 


judicious  manner.  It  turns  out  that  the  architecture  of  a  finite-element  program 
based  on  the  p-version  or  h-p-version  will  be  significantly  different  from  that  of 
one  of  the  standard  programs  founded  on  the  h-version.  Moreover,  the  mathematical 
properties  of  these  different  versions  are  rather  distinct  as  well. 

In  order  to  illustrate  the  last  point,  we  mention  a  typical  rate— of -convergence 
result  proved  in  [10] ,  [11] .  Suppose  that  the  use  of  the  h-version  with  a  sequence  of 
uniform  meshes  leads  to  a  rate  of  convergence  of  o'der  N-01  where  a  does  not  involve 

the  (fixed)  degree  of  the  elements.  There  the  p-version,  in  general,  exhibits  a 

—6 

rate  of  order  S  where 

a  -  c  <Si2a-e 

with  any  e  >  0. 

For  the  h-version  the  convergence  rate  is  governed  by  the  degree  of  the  elements 
and  the  smoothness  of  the  solution.  More  specifically,  for  a  wide  class  of  practical 
problems  involving  only  singularities  caused  by  comers  of  the  domain,  etc.,  it 
can  be  proved  (see  [12])  that  there  exist  sequences  of  suitably  refined  (non-uniform) 
meshes  for  which  the  rate  of  the  h-version  is  controlled  solely  by  the  degree  of 
the  elements.  On  the  other  hand,  in  [11]  it  was  shown  that  in  the  same  cases,  the 
h-p-version  provides  for  a  convergence  rate  which  is  better  than  any  polynomial  in 
N  or  even  exponential. 

In  general,  the  p-version  is  better  in  handling  singularities  (located  at 
the  vertices  of  elements’!  than  the  h-version  does  in  the  case  of  uniform  meshes. 
Moreover,  it  has  been  shown  in  [13],  [14],  and  [15]  that  the  p-version  is  practi¬ 
cally  not  influenced  adversely  by  the  presence  of  nearlv  incompressible  materials. 

On  the  other  hand,  the  h-version,  especially  when  low  order  elements  are  used,  has 
encountered  difficulties  due  to  this  incompressibility.  In  many  cases,  the  use  of  the 
so-called  selectively  reduced  integration  techniques  helps  to  overcome  these  diffi¬ 
culties.  These  techniques  are  equivalent  to  some  mixed  method.  For  a  typical  analysis 
of  problems  of  this  type  in  one  dimension  we  refer  to  [16] . 
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These  comments  already  indicate  that  future  developments  should  take  more 
account  of  the  selection  of  the  best  form  of  the  finite  element  method  for  particular 
classes  of  computational  problems.  In  particular,  the  h-p-version  appears  to  offer 
an  attractive  opportunity.  It  is  certainly  beginning  to  be  recognized  now  that  the 
use  of  high-order  elements — when  properly  implemented  in  a  hierarchial  manner — is 
more  effective,  in  many  cases,  than  the  classical  h-version  approach.  For  some 
related  references  we  refer  to  f 1 3] ,  [17]  and  references  there.  The  p-version  as 
well  as  the  h-p-version  were  implemented  in  the  commercially  available  program 
COMET -X  for  2  and  3  dimensional  problems.  But  it  remains  a  research  problem  how 
to  design  effective  adaptive  approaches  for  the  h-p-version. 


/ 


3 .  A  Posteriori  Estimates  and  Adaptive  Techniques. 

The  standard  results  about  errors  of  finite  element  solutions  concern 
almost  exclusively  a  priori  estimates.  As  noted  before,  as  important  as 
these  estimates  are  for  theory  cf  the  method,  they  are  rarely  useful  in 
practice  for  estimating  the  error  of  a  particular  computed  solution.  For  that 
we  have  to  turn  to  a  posteriori  estimates  which  utilize  information  obtained 
during  the  solution  process  itself. 

The  use  of  a  posteriori  estimates  means  that,  together  with  the  approximate 
solution,  the  error  estimator  c  is  computed  which  represents  a  measure  of  norm 
| j  e  j }  of  the  error  e  of  that  solution.  As  a  measure  of  the  reliability  of  the 
estimate  we  may  introduce  the  effectivity  index 


e  = 


Clearly  6  should  be  close  to  one  when  the  error  [ | e |  j  diminishes. 

We  distinguish  estimators  by  the  information  about  0  that  can  be 
guaranteed.  An  error  estimator  c  is  an  upper  estimator,  or  a  lower  estimator 
if 

llell£  ^ 

or 

CL  e  <  I |e| I  , 

respectively,  where  the  constants  depend  on  reasonably  large  class  of  possible 
solutions  and  on  a  (reasonably  large)  class  of  admissible  finite  element  meshes, 
but  not  on  the  specific  computed  solution  and  its  mesh.  An  upper  or  lower 
estimator  is  called  "guaranteed"  if  C^=l  or  0^=1  ,  respectively.  Finally,  an 
estimator  is  asymptotically  correct  if 


lim 


0 


c 


1 


for  a  specified  class  of  admissible  meshes.  In  general,  this  class  turns  out 
to  be  more  restrictive  than  that  used  in  any  property  of  c  to  be  an  upper  or 
lower  estimator.  For  an  analysis  of  these  properties  of  various  estimators 
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for  a  class  of  problems  in  one  space  dimension  and  for  different  norms,  we 
refer  to  [ 18] . 

For  practical  application  it  is  essential  that  1 1*—  1 )  is  small,  say  of  the 
order  of  0.2  to  0.3  when  the  relative  error  is  of  the  order  of  10%  and  the 
order  of  0.05  to  0.15  when  it  is  in  the  range  of  2-5%.  In  general,  the  size 
of  i 6 — 1 [  is  a  much  more  important  item  of  information  than,  say,  the  fact  that 
C^.=l  ,  that  is,  that  ,  is  a  guaranteed  upper  estimator.  Theoretical  results 
as  well  as  experimental  experience  indicate  that  the  type  and  method  of  construc¬ 
tion  of  the  finite  element  meshes  is  crucial  for  | 6—1 ]  to  be  small.  Not 
surprisingly,  these  classes  of  meshes  can  only  be  obtained  by  adaptive  techniques. 

The  classical  finite  element  method  addresses  positive  definite,  self- 
adjoint  problems.  Then  the  energy  norm  is  the  natural  norm  even  though  it  is 
not  always  the  most  relevant  one  for  practical  applications.  When  classical 
forms  of  the  finite  element  method  are  considered,  the  theory  based  on  the 
energy  no: m  is  especially  advantageous  and,  for  instance,  the  well-known  two- 
sided  estimates  of  the  error  becomes  available.  But  their  practical  computation 
is  expensive,  and,  as  all  such  global  constructions  of  error  estimators,  they 
turn  out  to  be  unusable  as  a  tool  for  adaptive  mesh-constructions.  In  other 
words,  we  should  look  for  error  estimators  e  with  a  local  character.  By  this 
we  mean  that  t  is  composed  of  error  indicators  n(x)  each  of  which  corresponds 
to  a  single  element  x  of  the  current  mesh  A  and  is  computed  from  data  about 
the  solution  on  t  and  at  most  the  immediate  neighbors  of  this  element.  The 
construction  of  the  error-indicator  o(t)  and  their  composition  into  the  overall 
estimator  e  depends  on  the  particular  norm.  For  example,  for  the  energy  norm 
one  is  led  to  the  definition 

£  «  [  7.  n(A)2]'2  . 

t 

All  the  estimators  studied  in  [18]  and  [19]  for  the  one-dimensional  and  in  [20] 
for  the  two  space  dimensions  have  been  local  character. 
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The  local  nature  of  these  estimates  allows  their  use  in  the  construction 
of  meshes  with  optimal  error  behavior.  There  are  various  papers  in  the 
literature  (see  e.g.  [21],  [22],  [23])  which  address  the  characterization  of 
optimal  meshes  for  a  prescribed  number  of  nodal  points  and  a  specified  norm. 

From  a  theoretical  viewpoint  these  results  are  certainly  important,  since  they 
provide  information  about  the  best  possible  error.  But  for  practical  applica¬ 
tions  the  question  of  looking  for  an  optimal  mesh  is  not  very  effective.  In 
fact,  we  are  interested  only  in  achieving  an  error  that  is  nearly  optimal  and 
it  turns  out  that  there  is  considerable  flexibility  in  the  meshes  for  which 
this  holds. 

This  and  related  results  were  presented  and  analyzed  in  depth  in  [19] 
for  a  class  of  one-dimensional  problems.  There  it  also  showTn  that  asymptotically 
(in  a  certain  mathematical  sense  of  the  word)  the  meshes  with  optimal  error  are 
characterized  by  the  equality  of  the  error  indicators  p(t)  mentioned  above. 

This  represents  the  so-called  equilibration  principle  for  mesh-construction. 

It  also  turns  out  that  this  equilibration  is  essential  to  guarantee  that  6 
leads  to  one  (see  [20]). 

Although  the  results  about  the  equilibration  have  been  fully  proved  only  for 
one-dimensional  problems,  all  indications  are  that  in  some  form  or  other  they  will 
hold  in  general.  In  any  case,  we  should  not  expend  too  much  costly  effort  in  the 
determination  of  an  optimal  mesh,  but  instead  construct  a  succession  of  meshes 
which  are  as  much  as  possible  equilibrated  and,  hence,  provide  for  solutions 
w»ith  nearly  optimal  errors  for  the  corresponding  numbers  of  degrees  of  freedom. 

The  algorithms  for  the  adaptive  construction  of  the  desired  sequence  of 
nearly  equilibrated  meshes  presently  are  still  largely  of  a  heuristic  nature, 

(see  e.g.  [3],  [24],  [25]  for  the  beginnings  of  a  theory  of  such  algorithms.)  In 

principle,  at  each  stage  of  the  process  any  element  t  of  the  current  mesh  A 
is  refined  for  which  the  error  indicates  n(t)  exceeds  some  threshold  and  its 
refinement  is  predicted  to  be  within  the  expected  equilibrations  range  of  the 
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resulting  mesh.  Since  a  new  solution  has  to  be  computed  on  the  refined  mesh 
which  is  expensive,  there  is  a  need  for  balancing  tha  quality  of  the  meshes  with 
the  effort  of  designing  them.  In  fact,  practical  interest  centers  not  only  on 
near  optimal  errors  for  the  number  of  degrees  of  freedom  used,  but  for  the 
total  machine  time  expended.  This  corresponds  also  to  the  stopping  criterion 
which  should  terminate  the  refinement  process  if  either  the  required  accuracy 
has  been  achieved,  or  the  machine-time  and  storage-estimations  are  exceeded. 

This  is  another  reason  in  favor  of  working  with  sequences  of  nearly  equilibrated 
meshes,  since  then  at  each  stage  the  error  of  the  computed  solution  is  expected 
to  be  nearly  optimal  for  that  stage. 


11 


4 .  An  Adaptive  Finite  Element  System 

Although  the  mathematical  theory  must  he  the  basis  of  understanding  the  con¬ 
cepts  discussed  in  the  previous  sections,  it  is  essential  to  test  them  experi¬ 
mentally  on  a  meaningful  set  of  realistic  problems.  In  order  to  facilitate  such 
an  experimental  program  a  prototype  system  incorporating  these  ideas  was  developed 
at  the  University  of  Maryland  in  recent  years.  It  was  dubbed  FEARS,  short  for 
Finite  Element  Adaptive  Research  Solver.  Details  about  the  system  may  be  found 
in  [24],  [25],  [26],  [27]  and  [28]  . 

FEARS  was  never  intended  to  become  a  production  system.  Instead,  in  its 
design  we  were  guided  by  the  following  axioms: 

(1)  To  verify  that  the  general  concepts  outlined  earlier  can  be  effectively 
implemented  in  a  robust  system. 

(2)  To  study  the  influence  of  the  adaptive  environment  upon  finite-element 
program  architectures,  the  relevant  data  structures,  etc. 

(3)  To  assess  some  potentialities  for  procedural  parallelity  in  the  method. 

FEARS  is  designed  to  solve  a  class  of  elliptic  systems  of  two  partial  dif¬ 
ferential  equations  in  two  dimensions.  More  specifically,  with  the  notation 


T 

and  similarly  for  v  =  (v^  ,  v^)  ,  the  admissible  problems  are  assumed  to  be 

written  in  the  weak  form 


B(u,v) 


fpu.T  .  3v  .  ,9u.T  _  TCT  3v  .  T_  1 

A^+  Bv  +  u B  -^  +  u CvJ 


dx^dx,, 


+ 


Ta 

u  Aqv 


ds 


JJ(D^  +  Ev)dxldx2 

n 


+ 


T  , 
e  v  ds 


1 

V 
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Here  A,  B,  C,  are  piecewise  constant  matrices  while  the  matrices  D,  E  and 
the  vector  e  are  given  functions  of  the  independent  variables  ■  The 

matrices  are  supposed  to  be  such  that  the  problem  is  elliptic  and  the  finite 
element  procedure  leads  to  symmetric  linear  equations. 

_  9 

The  domain  ..  C  FT  is  assumed  to  be  the  union  of  a  relatively  small  number 
of  curvilinear  quadrilaterals,  called  2D  subdomains.  Here  we  were  influenced 
by  the  successful  substructuring  concept  in  engineering.  The  number  of  2D  sub- 
domains  is  limited  by  storage  consideration;  in  practice,  we  did  not  exceed  the 
number  20.  The  (relatively  open)  sides  of  the  quadrilaterals  are  called  1D- 
subdomains  and  the  vertices  OD-subdomains .  The  ID-subdomains  are  assumed  to  be 
circular  arcs,  each  characterized  by  its  vertices  and  curvature.  The  union  of 
all  2D-,  ID-  and  OD-subdomains  is  the  set  F  C  fi  used  in  the  definition  of  the  weak 
form  B(u,v)  . 

The  matrices  A,  B,  C  are  constant  on  each  2D-subdomain  while  A^  is 
defined  on  T  and  constant  on  each  lD-subdomain.  Natural  boundary  conditions 
may  be  introduced  via  the  vector  e  . 

Since  each  2D-subdomain  is  assumed  to  be  dif feormorphic  to  the  unit  square 

9 

in  R“  ,  the  finite  element  meshes  on  the  2D-subdomains  are  defined  first  on  the 
unit  square  and  then  mapped  over  into  the  corresponding  curvilinear  quadrilateral. 

A  typical  mesh  on  the  square  is  shown  in  Figure  1,  it  consists  of  conforming 
bilinear  elements  of  different  size.  The  vertices  marked  by  crosses  are  irregular 
nodal  points.  The  values  of  the  solution  in  these  points  are  determined  by  the 
conformity.  The  mapping  into  tne  quadrilateral  is  constructed  by  blending  tech¬ 
niques  from  the  formulas  of  the  circular  side-segments  of  that  quadrilateral. 
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Figure  1 

The  admissible  error  norms  in  FEARS  have  the  form 


where  Q  is  an  input  matrix  and  e  the  difference  between  the  exact  and  the 
finite  element  solution.  The  choice  of  Q  may  be  different oneach  2D-subdomain. 
For  the  energy  norm  in  absence  of  B  and  C  we  use  Q  =  A  . 

The  error  indicators  n(x)  are  computed  for  each  element  of  the  mesh  A 
of  ft  as  functions  of  the  following  quantities: 

(1)  The  size  of  the  element  x  , 

(2)  The  jump  of  the  derivatives  across  the  boundary  of  t  ; 

(3)  The  projection  of  the  residual  of  the  computed  solution  on  t  onto 
the  space  of  constant  functions  on  x  . 

As  mentioned  in  the  previous  section  the  error  estimate  e  of  the  solution  is 
then  composed  of  the  individual  n(x)  .  For  further  details  about  the  error 
estimation  procedure  and  its  theory  we  refer  to  [20] . 

In  order  to  accomodate  the  type  of  meshes  and  the  refinement  process  a  spe¬ 
cial  data  structure  was  designed  in  the  form  of  a  rooted  tree.  A  detailed 
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description  of  this  design  and  the  corresponding  access  algorithms  is  given  in 
[28],  Tne  data  structure  provides  for  a  natural  pivot  ordering  based  on  nested 
dissection  for  the  solution  of  the  global  stiffness  matrices  by  LU-decomposition. 

The  adaptive  mesh-refinement  algorithm  follows  the  general  outline  given  in 
the  previous  section.  An  element  in  the  base-mesh  on  the  unit  square  is  refined 
by  dividing  it  into  four  congruent  subsquares.  The  threshold  used  in  deciding 
which  elements  are  to  be  subdivided  is  chosen  to  be  the  predicted  largest  error 
indicator  value  in  the  next  mesh.  This  prediction  is  based  on  a  simple  extrapo¬ 
lation  of  the  change  of  the  error  indicators  from  the  previous  to  the  current 
mesh.  The  threshold  may  be  revised  upward  if  the  resulting  storage  or  machine¬ 
time  requirements  are  predicted  to  exceed  the  available  resources,  or  the  desired 
accuracy  is  expected  to  be  reached. 

In  order  to  cut  down  on  the  cost  of  solving  too  many  large  lineai  systems, 
FEARS  employs  two  different  types  of  refinement  steps,  called  long  and  short 
passes,  respectively.  A  long  pass  involves  a  recomputation  of  the  full  solution 
on  the  entire  mesh  while  in  the  short  passes  only  certain  local  computations  are 
carried  out  which  allow  for  a  decision  about  further  refinements.  The  short 
passes  avoid  the  necessity  of  computing  a  solution  or  a  mesh  which  differs  only 
in  a  few  elements  from  the  previous  one.  The  combination  of  short  and  long  passes 
is  controlled  again  adaptively. 
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5 .  Some  Sample  Computations  with  FEARS 

During  the  past  year  FEARS  has  been  applied  to  an  extensive  range  of  practi¬ 
cal  problems  and  test  problems.  In  the  selection  of  the  test-problems,  we  were 
guided  by  the  desire  to  establish  a  worthwhile  set  of  benchmark  problems.  Such 
a  set  should  include  not  only  problems  for  which  the  assumptions  of  the  under¬ 
lying,  available  mathematical  theories  are  valid,  but  also  some  where  they  fail. 

In  assessing  the  outcome  of  any  particular  computational  experiment,  we  used 
a  number  of  criteria,  including,  in  particular  the  following  two: 

(1)  Computation  or  estimation  of  the  effectivity  index  0  as  a  function 
of  the  relative  accuracy  of  the  solution  for  the  given  norm. 

(2)  Calculation  of  the  observed  rate  of  convergence  as  a  function  of  the 
number  of  degrees  of  freedom  and  comparison  with  the  theoretically 
expected  rate  for  the  particular  type  of  problem. 

It  is  impossible  here  to  give  more  than  a  few  typical  examples  of  our 
experimentation  with  FEARS.  Those  included  here  were  clearer  mainly  as  an 
illustration  of  the  capabilities  of  the  system  rather  than  in  support  of  any 
specific  point  about  the  concepts  discussed  here. 

Example  1:  Consider  the  planar  elasticity  problem  for  the  ring  shown  in  Figure  2 
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Figure  2 
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with  prescribed  displacements 


u  =  1 ,  v  =  0  on  f 


u  =  v  =  0 


on  its  boundary  lines.  The  material  is  assumed  to  be  isotropic  and  homogeneous 
with  modules  of  elasticity  E  =  1  and  variable  Poisson  ratio  v  .  The  exact 
solution  of  this  plane-strain  problem  is  known  analytically.  The  error  of  the 
approximate  solution  in  the  energy  norms  is  given  by 

I  *e^E  “  (WFE  "  WEX^ 

where  W  and  W  denote  the  energy  of  the  exact  and  the  finite  element 
EX  FE 

solution,  respectively. 

Since  the  problem  has  various  symmetries,  only  a  quarter  of  the  ring  needs 
to  be  used.  In  FEARS  this  domain  is  taken  as  one  subdomain.  The  following 
Table  1  exhibits  some  experimental  data  for  this  problem  and  the  following  cases 

r!  =  1  ,  r0  =  -5  ,  and  rQ  =  .1 
v  =  0.3  ,  v  =  0.49  ,  and  v  =  0.499 
The  column  headings  are  as  follows: 

(1)  The  inside  radius  rQ 

(2)  The  Poisson  ratio  v 

(3)  The  number  of  elements 


(4)  The  number  of  degrees  of  freedom 

(5)  The  relative  error  in  percent  of  the  exact  solution:  100  ||ej|  /||u. 


EX1  'E 


(6)  The  computed  relative  error  100  e/ j  | u 

(7)  The  effectivity  index  8 
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Figures  3a-3f  show  sequences  of  meshes  adaptively  constructed  by  FEARS  in 


these  cases. 


The  cases 


are  denoted  as  follows: 


Fig- 

ro 

V 

3a 

•5 

.  3 

3b 

.1 

.3 

3c 

.  5 

.49 

3d 

.i 

.49 

3e 

.  5 

.499 

3f 

.1 

.499 

The  table  shows  clearly  the  high  quality  of  the  error  estimators 
independently  of  the  value  of  v  .  Clearly  | 9  —1 j  satisfies  the  earlier 
stated  requirements.  The  table  also  indicates  the  increasing  inef fectivity  of 
the  h-version  when  v  -  0.5,  that  is,  where  the  material  becomes  more  and  more 
incompressible.  Figure  3  shows  that  the  adaptive  procedure  refines  the  meshes 
in  a  very  plausible  manner.  This  point  will  become  more  evident  in  the  next 
example.  The  meshes  also  show  certain  small  anomalies.  Due  to  the  threshold- 
algorithm  mentioned  in  the  previous  section,  elements  someplace  miss  barely  in 
being  refined. 

Exanple  2:  Consider  the  plane-strain  problem  of  a  cracked  panel  as  shown  in 
Figure  4.  The  material  is  again  assumed  to  be  isotropic  and  homogeneous  case 
E»1  and,  plane  strain,  and  Poisson  ratio  v  =  0.3  .  Again  because  of  symmetry 
we  may  work  only  with  the  upper  right  quadrant. 
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Figure  6 


N  *  101  N  «  180 

Figure  5  shows  the  accuracy  of  the  solution  in  terms  of  the  error  in  the 
energy  (not  the  energy  norm)  for  the  adaptively  constructed  meshes  given  in 
Figure  6.  In  addition,  Figure  5  includes  a  comparison  of  the  performance  of 
the  h-version  and  the  p-version.  More  specifically,  for  a  partition  of  the  domain 
into  eight  congruent  triangles  the  p-version  clearly  exhibits  twice  the 
rate  of  convergence  as  the  h-version  for  uniform  meshes.  This  agrees  with  the 
general  result  mentioned  in  Section  2.  We  also  observe  the  earlier  noted  best 
possible  rate  of  order  N  2  for  the  sequence  of  adaptively  constructei  meshes. 

Of  course,  the  corresponding  energy  rate  is  0(N  . 

In  problems  of  this  type  the  stress-intensity  factor  is  an  important 
quantity  for  engineering  considerations.  It  represents  a  partial  derivative 
with  respect  to  the  crack  length.  Accordingly,  an  error  estimator  of  the  stress 
intensity  factor  can  be  obtained  from  a  derivative  of  the  error  estimators  of  the 
solutions.  It  can  be  shown  that  under  suitable  conditions  and  for  the 
sequence  of  adaptively  constructed  meshes  these  error  estimators  of  the  stress 
intensity  factor  are  asymt  otically  correct. 


In  order  Co  test  this,  the  boundary  conditions  of  the  present  problem  were 
adjusted  by  a  proper  boundary  condition  on  the  sides  A^A?  and  A^A-,  so  that  the 
solution  is  known.  Table  2  shows  some  results  computed  with  FEARS  for  this  case. 


N 

Solutio 

rel.  t:rror 

n  error 

6 

Error  of  stress 

rel.  error 

-intensity  factor 

e 

22 

21.4% 

0.59 

12.0% 

0.38 

52 

12 . 6% 

0.74 

4.3% 

0.60 

109 

7.6% 

0.85 

1.7% 

0.76 

238 

4.7% 

0.93 

0.6% 

0.95 

Table  2. 


As  the  table  indicates,  the  rate  of  convergence  of  the  stress  intensity  factor 
is  twice  as  high  as  that  of  the  solution.  Under  appropriate  conditions  this 
can  be  proved  theoretically.  We  also  note  that  the  effectivity  index  of  the 
stress  intensity  factor  is  approximately  the  square  of  that  for  the  solution. 
Example  3:  As  a  final  example  we  consider  a  seepage  problem  involving  a 
piezometric  head  u  which  sa  isfies  the  equation 

2 

:  a  pL.  „  o 

.  .  O  x.  ex. 

1=1  l  i 

where  a  represent  the  permutabilty  and 


are  the  velocity  components.  Figure  7  shows  the  domain  partitioned  into 
sixteen  2D-subdomains .  In  each  one  of  these  subdomains  the  permeability  is 
assumed  to  be  constant.  The  values  of  the  constants  are  speficied  in  Table  3. 

As  is  typical  for  problems  of  this  kind,  we  set  u  *=  X2  at  the  upper  boundary 

c_u 
3  n 


and  use 


0  along  the  vertical  and  lower  sides  of  the  domain. 


E 


Figure  7 


sub 

domain 

a 

sub 

domain 

a 

1 

5.0 

9 

1.0 

2 

5.0 

10 

5.0 

3 

5.0 

11 

5.0 

4 

5.0 

12 

5.0 

5 

1.0 

13 

60.0 

6 

2.0 

14 

1000.0 

7 

70.0 

15 

5.0 

8 

1000.0 

16 

5.0 

Table  3. 
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For  the  error  assessment  the  following  two  norms  were  used 
(1)  The  energy  norm 

/  [  I  ft  v’  J.  1 

U 


J 


1  a  (i - )  dx.  d>  „ 

.  ,  dx  .  1  £ 

1=1  1 


(2)  The  velocity  norm 
2 


IMtv  =  (  |  [*  vi2  dxl  dx2  ) 
•  ' . -  1=1  / 


Table  4  gives  some  experimental  data  obtained  with  FEARS  for  this  problem. 
(The  effectivity  index  0  is  only  available  for  the  energy  norm).  (The  exact 
energy  was  computed  by  still  larger  computations  and  various  extrapolation 
techniques,  so  that  the  true  error  can  be  established.) 


No.  of 

mi 

Norm 

elements 

mSEIBIR 

0 

1  1  •  I  1 

451 

8.16 

1.012 

II  1  lE 

850 

6.7706 

6.92 

0.982 

1  1  .  1  1 

478 

6.8035 

N.  A 

II  1  1  v 

886 

6. 7776 

N.  A. 

Table  4. 


The  adaptive  approach  based  on  the  energy  norm  is  expected  to  perform 
better  than  that  based  on  the  velocity  norm  when  energy  norm  of  the  errors  is 
considered.  This  agrees  with  the  data  in  Table  4.  The  corresponding  number  of 


elements  in  the  sixteen  subdivisions  is  given  in  Table  5. 
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sub 

domain 

Energy 

— 

norm 

- - - 

Velocity 

norm 

1 

19 

28 

40 

46 

2 

55 

100 

55 

79 

3 

16 

37 

16 

28 

4 

10 

19 

34 

34 

5 

10 

25 

4 

22 

6 

10 

34 

4 

70 

7 

4 

4 

7 

25 

8 

7 

7 

31 

52 

9 

31 

97 

22 

25 

10 

76 

115 

55 

100 

11 

55 

91 

40 

58 

12 

64 

103 

40 

130 

13 

13 

19 

31 

55 

14 

4 

19 

31 

73 

15 

19 

28 

16 

19 

16 

58 

124 

52 

70 

Total  N 

451 

850 

478 

886 

Table  5 


We  end  this  section  with  a  few  comments  about  the  overall  performance 
of  FEARS.  The  system  includes  a  set  of  timers  and  counters  to  measure  the 
cost  of  various  program-modules.  An  an  example,  Table  6  gives  some  relative 
lines  for  the  seepage  problem  of  example  3  in  the  case  of  the  mesh  with  830 
constructed  adapted  adaptivity  using  the  energv  norm. 


Module 


Refinement  and  related  tasks  2.74% 

Elimination  of  the  blocks  corresponding  to  the 

2D  subdomains  and  overall  back  substitution  51.81% 

Elimination  of  the  .blocks  corresponding  to  the 

lD-and  0D-  subdomains  26.14% 

Computation  of  the  error  indicators  19.31% 


Total  100.00% 


Table  6 
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The  separation  of  the  elimination  process  into  a  two-dimensional  and  lower¬ 
dimensional  level  corresponds  to  the  decompostion  algorithm  based  on  the  earlier 
mentioned  tree-structure  for  the  data. 

The  table  provides,  of  course,  only  some  information  about  the  relative 
machine  time  for  the  final  mesh  of  the  process.  Experience  has  shown  that  the 
solution  on  the  final  mesh  may  be  expected  to  take  about  50%  of  the  total 
machine  line  for  the  overall  procedure.  Moreover,  in  our  runs  with  FEARS  we 
found  that  the  time  of  the  solution  on  the  final  nesh  is  order  n\  1.5  X  1.9, 
where  N  is  the  final  number  of  degrees  of  freedom.  If  one  takes  into  account 
that  the  final  mesh  provides  a  considerably  higher  accuracy  than  any  a-priori 
constructed  mesh  with  the  same  number  of  degrees  of  freedom  N  ,  then  the  overall 
efficiency  of  the  process  becomes  readily  apparent. 
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6 .  A  Convectlon-Dif fusion  Problem 

So  far  essentially  all  our  discussions  concerned  self-adjoint  linear  prob¬ 
lems.  For  the  remainder  of  this  survey  we  turn  now  to  problems  of  a  different 
nature. 

Convection-diffusion  problems  are  important  in  many  applications.  In  one 
dimension  a  typical  form  of  such  problems  is  as  follows: 

eu"  +  a(x)  u'  +  b(x)  u  =  f(x)  ,  0  <  x  <  1 

u(0)  =  u  ,  u(l)  =  6 

It  is  well-known  that  the  solution  shows  boundary  layer  effects  when  e  is  small 
either  because  of  the  boundary  conditions,  non-smoothness  of  f  ,  or  changing 
sign  of  a  .  Hence  it  is  not  surprising  that  the  numerical  solution  of  such 
problems  is  not  easy  and  runs  into  typical  difficulties  for  small  e  .  In  the 
design  of  solution  procedures  a  major  concern  is  its  robustness  and  stability 
with  respect  to  e  .that  is,  a  uniformity  of  performance  for  all  t  in  the  in¬ 
terval  0  <  e  <  1  .  In  addition,  in  the  application  of  the  finite  element  method 
special  care  must  be  taken  in  the  selection  of  the  trial  and  test  functions. 

Because  of  their  boundary  layer  character  problems  of  this  type  use  some 
form  of  the  norm  which  is  the  best  in  various  applications.  More  specific¬ 

ally,  in  connection  with  the  finite  element  method  it  becomes  necessary,  for 
theoretical  reasons,  to  introduce  the  following  mesh  dependent  L1  -norm: 


where  Xq  =  0  <  x^  <  . .  .  <  xN  <  xN+^  =  1  are  the  nodal  points  and  tr  « 

Xj  -  >  j  “  1,...,N+1  ,  the  mesh-steps.  For  details  about  these  norms  we 

refer  to  [29],  [30].  In  [31]  it  was  shown  for  the  case  of  piecewise  linear 
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test  functions  that  the  desired  uniformity  of  performance  with  respect  to  e 
requires  the  use  of  trial  functions  which  depend  on  e  ,  the  values  of  the  coef¬ 
ficient  functions  a, b  in  the  elements  and  their  sizes. 

The  use  of  these  functions  can  be  combined  with  a  posteriori  error  estima¬ 
tions  under  f  -norms  and  the  corresponding  adaptive  approaches.  Once  more 

it  can  be  shown  that  these  estimators  are  asymptotically  correct.  For  details 
we  refer  to  [32]  and  further  forthcoming  papers. 

As  an  illustration  consider  the  problem 

-eu"  +  v(x)  u' (x)  *  0  ,  0  <  x  <  1 

u(0)  =  2tah(-|)  ,  u(l)  =  0  , 


with 


v(x)  =  2 1 ah  ,  0  ^  x  <_  1 

This  is  a  model  of  Burger's  equation 

-eu"  +  uu'  =  0 

with  the  exact  solution  u  =  v  .  Figure  8  shows  the  accuracy  of  the  computed 


Figure  8 


solutions  for  uniform  meshes  as  well  as  adaptively  constructed  meshes  for  varii  us 


values  of  t  .  In  all  cases  the  L  -norm  was  used.  The  drastic  difference 

1  ,  u 

between  the  performance  in  the  two  cases  is  certainly  very  evident. 


Figure  9  presents  the  so-called  grading  function  S(x)  of  the  final,  adap¬ 


tively  constructed  mesh  (see  eq .  [18],  [32]). 


Figure  9 

The  relation  between  this  mesh  and  the  function  £  is  given  by 


^i>  "  I 

where  {x^}  again  denotes  the  nodal  points.  The  refinement  in  the  neighborhood 
of  x  =  1  where  the  boundary  layer  occurs  is  clearly  visible.  The  quality  of 
the  adaptively  constructed  mesh  sequence  is  apparent  from  Figure  9  where  the 
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observed  rate  of  convergence  is  0(N  ~)  which  is  the  maximally  possible  rate  for 
linear  elements. 

Finally  Figure  10  exhibits  the  behavior  of  the  effectivity  index  of  :he  es¬ 
timators  for  both  the  sequence  of  un. form  meshes  and  the  adaptively  constructed 
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Figure  10 


Again  the  high  quality  of  the  estimators  as  well  as  that  of  the  adaptive  mesh 
sequence  is  certainly  evident.  But  even  for  the  uniform  meshes  the  estimators  are 
not  too  bad.  Moreover,  the  figure  illustrates  the  theoretically  supported  fact 
that  the  estimator  is  simultaneously  an  upper  and  a  lower  one. 

The  general  layout  used  for  the  program  in  this  case  is  similar  to  that  of 
FEARS;  in  particular  the  adaptive  procedure  corresponds  entirely  to  that  employed 
here. 

Although  the  model  problem  is  only  one-dimensional,  the  ideas  formulated 


here  are  not  restricted  to  that  case  and  we  expe:t  that  they  can  be  applied  also 
to  more  general  situations  in  higher  dimensions. 
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7 .  A  Parabolic  Problem 

Among  the  most  effective  methods  for  solving  parabolic  problems,  as,  for 
instance,  che  transient  diffusion  problem 


is  the  method  of  lines.  Its  basic  idea  is  to  discretize  the  equation  in  the 
space  variables  and  to  solve  the  initial  value  problem  for  the  resulting  systeii 
of  ordinary  differential  equations  by  one  of  today's  sophisticated  adaptive 
ODE-solvers. 

If  the  finite  element  method  is  used  for  the  space  discretization,  then 
the  concepts  discussed  in  che  earlier  sections  suggest  the  following  tasks: 

(1)  Develop  strategies  for  the  adjustment  of  the  tolerance  in  the  ODE- 
solver  so  that  the  errors  caused  by  the  discretizations  in  space  and 
time  are  comparable. 

(2)  Design  and  analyze  error  estimators  for  the  total  error. 

(3)  Develop  algorithms  for  combined  space  and  time  discretization. 

Two  approaches  to  the  third  task  have  been  considered.  The  first  one  consists 
in  moving  the  nodal  points  of  the  finite  element  mesh  in  line  with  information 
obtained  during  the  solution  of  the  ordinary  differential  equations.  This  was 
discussed  in  [33],  [34],  [35].  In  the  second  approach  the  mesh  remains  fixed  during 
certain  intervals  in  time.  The  length  of  these  intervals  is  chosen  adaptively 
and  at  the  transition  between  them  the  mesh  is  changed  discontinuously ,  for  in¬ 
stance,  by  refining  (or  de-ref ining)  it  as  in  FEARS.  For  further  details  we  re¬ 
fer  to  [36],  [37]. 

The  adaptive  control,  of  course,  requires  the  availability  of  an  appropriate 
theory  of  a  posteriori  error  estimates.  Such  a  theory  was  introduced  in  [35], 

[37],  [38].  It  leads  to  estimators  which  turn  out  to  be  upper  and  lower  estima¬ 
tors  that  are  asymptotically  correct  under  suitable  conditions. 


A  general  procedure  for  a  class  of  parabolic  equations  in  one  space  variable 
has  been  implemented  based  on  the  mentioned  error  estimators  and  the  second  of 
the  above  approaches  for  adapting  the  finite  element  mesh.  As  an  illustration 
we  present  some  computational  results  for  the  problem 

I?  =  £  a(x’  S  "  b(x)u  +  f(t’x)  ’  0  <  x  <  1  *  C  >  0 

u(t,x)  =  g(t,x)  ,  X  =  0,1  ,  t  >  0 
u(0,x)  =  Uq(x)  ,  0  <  x  <  1  • 

More  specifically,  let 

a(x)  =  cosh(4x-2)  ,  b(x)  =  sinh(2x) 


and  choose  uQ  ,  f  ,  and  g  such  that  the  exact  solution  becomes 


u(t,x) 


lOt  +  l/10y 


+  l-i(l+tan(10t-ll)  )  (2sin(nx)  +  2sin(2TTt)sin(2irx)) 


This  oscillatory  solution  was  selected  in  order  to  show  the  effectivity  of  the 
approach.  Figure  11  illustrates  the  form  of  the  solution  for  various  t  . 


X 


Figure  11 

For  the  computation  uniform  meshes  with  N  intervals  and  piecewise  linear 
elements  were  used  and  the  ODE-solver  LSODI  was  applied.  Figure  12  shows  the 
relative  errors,  measured  i:i  the  energy  norm 


N  •  10 


H  -  tO 


N  •  *0 


Figure  12 
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and  Figure  13  the  corresponding  effectivity  indices  as  functions  of  t  and  N  . 
Once  again  we  see  the  high  quality  of  the  estimators  even  for  relatively  modest 
accuracies . 


Figure  13 

The  interaction  between  the  device  of  tolerance  for  the  ODE-solver  and  the 
space  discretization,  as  well  as,  the  design  of  a  total  error  estimator  is  dis¬ 
cussed  in  [38].  The  results  can  be  extended  to  problems  in  higher  space  dimen¬ 
sion  and  also  to  mildly  nonlinear  parabolic  equations. 
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8.  Nonlinear  Problems 

In  recent  years  the  application  of  finite  element  techniques  to  the  numeri¬ 
cal  solution  of  nonlinear  problems  has  been  receiving  increasing  attention.  In 
particular,  the  nonlinear  analysis  of  the  behavior  of  stationary  mechanical  sys¬ 
tems  has  become  more  and  more  important  in  engineering  applications.  We  refer, 
for  example,  see  [1].  The  area  is  rather  large  and  we  can  present  here  only  a 
brief  overview  of  some  questions  and  results  relating  to  the  general  concepts 
discussed  earlier. 

8.1  Parametrical  Systems:  When  we  speak  of  the  behavior  of  a  mechanical  struc¬ 
ture,  then  we  mean,  for  instance,  the  deformation  of  the  structure  in  response 
to  the  action  of  external  influence  quantities  such  as  load  distributions,  changes 
in  material  properties  and  geometrical  data,  etc.  This  implies  that  the  equations 
which  govern  the  deformation  always  involve  a  set  of  influence  parameters,  and, 
in  line  with  the  nonlinear  character  of  these  equations,  we  cannot,  in  general, 
draw  conclusions  about  solutions  for  arbitrary  parameter  values  from  one  such 
solution. 

Hence,  after  suitable  finite  element  discretizations  we  are  led  to  consider 
finite-dimensional  nonlinear  equations  of  the  form 


F(y,p)  “  0 

where  y  6  Rn  is  a  vector  of  state  variables  (e.g.,  deformations),  p  €  R™  a 
vector  of  parameters,  and  F  a  given  function  which  maps  Rn  x  Rm  into  Rn  . 
The  study  of  the  response  of  the  corresponding  mechanical  system  now  requires  a 
determination  of  parts  of  the  solution  set 


E(F)  -  { (y ,p)  €  Rn  x  F(y,p)  =  0} 
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and-  a  study  of  its  shape  and  features.  In  structural  analysis  E(F)  is  often 
called  the  equilibrium  surface  of  F  . 

Suppose  now  that  F  is  of  class  Cr  ,  r  1  ,  and  that  the  derivative  of 

F  at  (y ,p)  £  Rn  x  Rm  is  denoted  by  DF(y,p)  .  If  the  regularity  set 

R(F)  =  {(y,p)  £  Rn  x  Rm;  rank  DF(y,p)  =  n} 

is  non-empty,  then  the  regular  part, 

Er(F)  =  E(F)  0  R(F)  , 

of  the  equilibrium  surface  of  F  is  an  open,  m- dimensional  Cr-manifold  in 

Rn  x  Rm  .  At  a  singular  point  (y,p)  €  E(F)  —  E  (F)  the  surface  may  show  a  vari- 

ety  of  different  features  but  we  will  not  discuss  this  here. 

The  characteristics  of  the  equilibrium  surface  depend  strongly  on  the  selection 
of  the  parameters.  Improperly  selected  parameters  may  well  introduce  meaningless 
singularities.  Thus  the  selection  of  physically  meaningful  parameters,  especially  in 
neighborhoods  of  singular  points  is  of  considerable  importance. 

From  a  computational  point  of  view  it  is  not  unexpected  that  the  larger  the 
set  of  free  parameters  the  more  involved  is  the  task  of  analyzing  the  equilibrium 
surface.  Hence,  in  practice,  the  standard  approach  is  to  restrict  the  variation 
of  the  parameter  vector  p  to  one  degree  of  freedom.  In  other  words,  we  require 
that  p  is  a  function  p  =  p(t)  of  one  real  variable  t  .  Then  the  given  equa¬ 
tion  reduces  to  an  equation 


class  Cr  ,  then  the  regular  part  E^(F)  of  the  equilibrium  surface  of  F  con- 
stitutes  a  one-dimensional  C  -manifold,  that  is,  a  curve  on  E(F)  . 

It  can  be  seen  that  we  have  (y,t)  £  R(F)  exactly  if  one  of  the  following 
two  conditions  holds: 

(i)  D^F(y,p(t))  is  non-singular. 

(ii)  rank  D^F(y,p(t))  =  n-1  and  D^F(y ,p ( t) )p ’ (t )  ^  rge(DyF(y,p(t)) 

In  either  case  we  have  (y,r(t))  £  R(F)  .  The  points  (y,t)  £  E„(F)  where  (ii) 
holds  are  called  limit  points  or  turning  points  of  F  with  respect  to  t  . 

On  its  regularity  set  R(F)  the  function  F  defines  a  direction  field. 

More  specifically,  there  exists  a  unique  mapping  T  from  R(F)  into  R0"*^ 
such  that 

/i>F<y,o\ 

DF(y.t)T(y,t)  «  0  ,  ||T(y,t)||,  -  1  ,  det  _  >  0  , 

\T(y,t)T/ 

(see  eq.  [39]).  Then  (y,t)  is  a  limit  point  of  F  with  respect  to  t  exactly 
if  the  (n+l)st  component  of  T(y,t)  is  zero. 

In  the  case  of  structural  problems,  the  critical  points  encountered  on  the 
solution  curve  defined  by  F  ,  that  is,  the  limit  points  of  F  or  the  points 
not  in  R(F)  ,  signify  a  possible  loss  of  stability  of  the  structure.  Such  a 
loss  of  stability  actually  corresponds  to  a  dynamic  phenomenon  whereby  the  struc¬ 
ture  undergoes  a  sudden  change  of  deformation.  The  dynamics  of  this  phenomenon 
are  not  specified  by  the  equations  involving  F  or  F  ,  but  it  is  possible  to 
deduce  from  the  shape  of  the  equilibrium  surface  what  type  of  sudden  changes  may 
be  expected.  Broadly  speaking,  at  limit  points  we  expect  to  encounter  a  snap- 
through  or  collapse  behavior  while  at  points  not  in  R(F)  the  structure  say 
undergo  (stable  or  unstable)  buckling. 
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These  observations  already  indicate  that  for  a  deeper  understanding  of  the 
behavior  of  a  structure  it  is  necessary  to  analyze  the  form  and  principal  features 
of  the  equilibrium  surface  of  F  .  For  this  we  need  to  consider  continuation 
methods  in  a  broader  sense  as  a  collection  of  numerical  procedures  for  accom¬ 
plishing  tasks  of  the  following  type 

(i)  Follow  numerically  any  curve  on  the  surface  specified  by  a  given  con- 
bination  p  =  p(t)  of  parameter  values. 

(ii)  Determine  on  such  a  curve  the  critical  points  when  stability  may  be  lost, 

(iii)  From  any  critical  point  follow  a  path  in  the  boundary  of  the  set  of 

stable  solutions  on  the  surface. 

(iv)  On  any  curve  specified  by  (i)  determine  the  location  of  points  not  in 
R(F)  and,  in  particular,  of  bifurcation  points  from  which  other  paths 
emanate,  and  trace  any  such  path. 

Methods  of  these  and  related  tasks  have  been  developed  by  many  authors  in  recent 
years  (see  eq.  [40],  [41],  [42],  [43],  [44],  [45],  [46],  [47], [48]).  We  shall 
not  go  into  details  here  but  sketch  in  the  next  section  only  some  of  the  principal 
features  of  methods  for  the  first  two  tasks  (i)  and  (ii) . 

The  discussion  in  this  section  certainly  shows  that  in  the  nonlinear  case 

the  problem  formulation  as  well  as  the  computational  requirements  differ  strongly 
from  those  in  the  linear  case.  In  addition,  there  is  now  not  only  a  need  for  es¬ 
timates  of  the  error  of  the  computed  solution  but  also  of  the  various  computed 
features  of  the  equilibrium  surface,  as  for  instance,  of  the  location  of  limit 
points  or  bifurcation  points,  etc.  Furthermore,  there  is  a  requirement  for  a 
combined  adaptive  control  of  the  space  discretization  and  the  steps  of  the  con¬ 
tinuation  process.  In  a  sense  this  is  similar  in  concept  to  the  situation  dis¬ 
cussed  earlier  in  the  case  of  parabolic  equations. 


As  we  saw  our  interest  centers  on  computing  continuous  paths 


x  =  x (s )  =  (y  (s) ,  t  (s) )  €  Rn+\  s  €  1 


on  E^(F)  where  F  is  a  reduced  mapping  of  the  form  discussed  in  the  previous 
section,  and  s  is  a  suitable,  as  yet  unspecified  parameter  ranging  over  some 
interval  I  €  .  For  ease  of  notation  we  denote  the  vectors  (y,t)  €  Rn  x 

simply  by  x  . 

For  any  x  €  R(F)  the  direction  vector  Tx  defined  earlier  may  be  computed 
as  follows 


DF(x) 

(eV 


v  =  e 


n+l 


sign  det\ 


Tx  =  o 


■ 


DF(x) 
(e  )  / 


i  T 

sign(e  )  v 


where  e1 , . . . , en+1  are  the  natural  basis  vectors  of  Rn+^  and  e1  was  chosen 

k  A 

such  that  Tx  has  a  component  in  its  direction.  If  x  €  R(F)  is  a  known 
approximation  of  our  derived  path  x  =  x(s)  then 


Tt(h)  =  x^  +  h  Tx^  ,  h  €  R^ 


k 

represents  an  approximation  of  the  tangent  line  of  the  curve  at  x  .  Hence,  if 
i  i  t  k 

e  again  is  such  that  (e  )  Tx  #0  then,  for  sufficiently  small  h  ,  the  hyper- 
i  T 

planes  (e  )  (x— rr(h))  =  0  are  transversial  to  the  line  tt  =  r(h)  and  hence  also 
to  the  curve  x  =  x(s)  itself.  In  other  words,  the  augmented  equation 


Sh(x) 


\ 

v(ei)T(x-7i(h  )y 


0 


has  a  unique  solution  for  any  fixed,  small  h  which  may  be  computed,  say,  by 
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Newton's  method  applied  to  this  equation  with  Ti(h)  as  starting  point. 

With  these  observations  we  may  outline  a  so-called  locally  parametrized 
continuation  process  as  follows: 

1.  Given  x°  £  E  (F)  and  an  initial  index  i  .  ,  set  k  =  0; 

K  —  i. 

k 

2.  Compute  Tx  using  the  index  i  =  i^  ^  ; 

i  T  V 

3.  Che  ose  i  =  i^  ,  1  <  <  n+1 ,  such  that  (e  )  TxK^  0; 

4.  Choose  a  steplength  h  along  the  line  r  =z(h); 

5.  With  w(h)  as  starting  point  apply  Newton's  method  to  g^(x)  =  0; 

k+1 

6.  If  "convergence"  then  x  =  computed  solution,  k  =  k+1,  go  to  2; 
else  go  to  5  with  h  =  %h  . 

In  line  with  this  schematic  outline  complete  programs  have  been  developed. 

For  a  specific  version  ve  refer  to  [49]. 

Because  of  the  significance  of  limit  points  in  many  applications,  it  is  not 

surprising  that  a  number  of  algorithms  have  been  proposed  for  the  computation 

of  such  points.  We  mention  here  the  methods  presented  in  [49],  [50],  [51],  [52], 

[53],  [54],  [55],  [56],  [57]  .  There  are  two  basic  classes  of  algorithms 

for  computing  limit  points,  namely  (a)  local  iterative  processes,  and  (b)  methods 

utilizing  the  continuation  path.  The  members  of  the  first  class  may  be  started 

* 

anywhere  in  the  neighborhood  of  the  desired  limit  point  x  while  those  in  the 
second  class  require  an  input  some  point  (or  points)  on  the  continuation  path 

ft 

through  x  .  In  practice,  this  distinction  tends  to  loose  its  significance 

since  also  in  the  first  case  a  continuation  process  is  needed  to  reach  a  neighbor- 
* 

hood  of  x 

In  [58]  experimental  comparisons  of  all  the  mentioned  methods  for  a  set  of 

test  problems  were  presented  and  some  guidelines  for  their  applicability  given. 

A 

It  turns  out  that  when  two  approximations  bracketing  the  limit  point  x  are 

'k 

available  on  the  continuation  path  through  x  then  a  version  of  the  method 
given  in  [49],  [54]  is  highly  reliable  and  efficient  and  hence  can  be  recommended. 
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8 . 3  A  Model  Problem 

Not  unexpectedly,  for  nonlinear  problems  the  theory  of  a  posteriori  error 
estimators  has  not  yet  advanced  as  far  as  in  the  linear  case.  However,  a  number 
of  results  have  been  developed  and  proved  already  and  for  a  detailed  survey  we 
refer  to  [39].  Here  we  cannot  enter  into  details,  but  restrict  the  discussion 
to  some  computational  results  for  a  specific  model  problem. 

The  problem  under  consideration  concerns  the  planar  deformations  of  a  non- 
linearly  elastic  rod.  More  specif ically,  large  deformations  as  well  as  material 
non-linearities  are  allowed.  A  comprehensive  treatment  of  the  global  behavior 
of  such  rods  was  given  in  [59].  These  rods  can  suffer  not  only  flexure  but  also 
compression  (or  extension)  and  shear.  The  basic  conf iguration  of  the  rod  is 
sketched  in  Figure  14. 


Figure  14 

We  restrict  ourselves  to  the  case  of  a  symmetric  shape  where  it  suffices  to  con¬ 
sider  only  one  half  of  the  rod,  say  the  left  half  in  the  figure.  Then  the  con¬ 
figuration  of  the  rod  is  described  by  two  functions 

r:  [0,1]  -  R2  ,  9:  [0,1]  -  R2  ,  ||q(s)||2  -  1 
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where  r  defines  the  axis  of  the  deformed  rod  and  q  the  direction  of  the  cross- 
1  0 

section.  Let  e  ,e“  be  the  global  basis  vectors  as  shown  in  the  figure  and  u 
the  angle  between  q  and  e~  .  Under  the  influence  of  a  constant  load  with  the 
components  X  and  .  as  indicated  in  the  figure,  the  problem  can  then  be  modeled 
by  the  two-point  boundary-value  problem 

-p-  a(u',u,X,v)  +  (Xsin  u  +  v  cos  u) (1  +  b(-Xcos  u  +  v  sin  a)) 
ds 

+  (Xcos  u  -  v  sin  u)c(Xsin  u  +  v  cos  u)=0,  0  <  s  <  1  , 

u(0)  =  u(l)  =  0  , 


with  given  functions  a,  b,  c  representing  constitutive  equations. 

For  our  computations  we  employed  the  standard  finite  element  method  with 
C°  elements  of  degree  p  . 

For  the  specific  choice  of  the  coefficient  functions 


a(t)  =  c(t)  =  arctan  ^t  ,  b(t) 


Figure  15  shows  several  solution  curves  for 
system  with  X  as  abscissa  and  the  norm  j  |u 


t 

1-t 

for 

t  <  0 

2 

t+2t 

for 

t  >_  0 

v  =  constant  plotted  in  a  coordinate 
| |  of  the  computed  solution  as  or- 


Loo 

dinate.  In  all  cases  uniform  meshes  with  eight  sub-intervals  and  cubic  C°-elements 
were  used. 

In  order  to  assess  the  effectivity  of  the  error  estimation,  we  used  a  "nearly 
exact"  solution  the  computed  solution  with  a  uniform  mesh  of  10  subintervals  and 
quintic  C°-elements.  We  consider  the  curve  v  =  0.1  and  on  it  the  five  target 
points  X  =  1 . 0, 2 . 5 , 3. 5 , 2 . 5 , 1 . 6 .  Table  7  shows  the  computed  error  estimators  of 
the  error  in  the  norm  and  the  H^-norm  of  the  "nearly  exact"  error  in  percent 
of  the  norm  of  the  solution.  The  table  shows  clearly  the  high  quality  of  the  es¬ 
timators.  In  the  cases  left  blank  in  the  table  the  error  is  essentially  zero 


within  machine  accuracy. 
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2  3 


1 

Error  estim. 

; 

12.15% 

0.1137% 

0.004089% 

1  .0 

Error  norm  ; 

i 

■  12.13% 

0.1138% 

0.004089% 

Error  estim.  j 

11 .55% 

0.4306% 

0.01552% 

1 

I 

2.5 

: 

i 

Error  norm 

12.04% 

0.4310% 

j  0.01552% 

1 

Error  estim. 

3.5 

j  21.58% 

1 .597% 

1 .123% 

0.5011% 

Error  norm 

-  - 

23.91% 

i 

1 . 900% 

1  .391% 

0.5389% 

Error  estim. 

2.5 

21 . 53% 

** 

CO 

o 

r-. 

1 .003% 

0.5013% 

Error  norm 

21 .88% 

-  -  - 

1.865% 

1 .285% 

0.5455% 

Error  estim. 

1 .6 

19.45* 

0.7632; 

0.8123.%  -- 

0.3171% 

Error  norm 

19.14% 

1.135% 

0.9576% 

0.331 2% 

As  Figure  15  shows  there  are  two  limit  points  along  the  curve  for  v  =  0.1. 

With  a  mesh  of  eight  equal,  cubic  C°-elements  these  limit  points  occur  at 

A  =  3.6546  and  •  =  1.46749.  Table  7  shows  the  effectivitv  of  our  error- 

1  *2 

estimators  of  these  critical  values.  For  the  comparison  five  and  ten  equal, 
quadratic  C°-elements  were  used.  Again  the  estimators  perform  very  well. 

In  general  the  error  varies  considerably  along  any  one  solution  curve. 


1 

!  £st.  of 

! 

Limit  Point 

L  i  mi  t 

No.  of  ! 

Solution 

| 

Point 

Interv. 

Error 

Value 

Error 

Estim. 

— - ■■■—"■  —  -r 

5 

13. ns 

3.699385 

0.04472 

0.05294 

\ 

10 

3.18% 

j 

j  3.657982 

i 

0.00317 

0.00399 

• 

5 

7.205% 

!  1.470290 

0.002793 

0.002843 

Xj2 

10 

1.613% 

1  .467897 

0.00040 

0.000357 

Table  8 

Thus  adaptive  approaches  are  extremely  important  in  order  to  control  this  error 
along  the  curve.  Our  experimentation  has  shown  that  savings  of  about  70%  or 
more  in  computational  effort  can  be  achieved  for  problems  of  the  type  discussed 
here.  For  higher  dimensional  problems  such  savings  become  extremely  valuable. 
For  more  details  we  refer  to  [39]. 
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